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, I Abstract. We revisit the problem of extending quadrature formulas for general weight 

C^ ■ functions, and provide a generalization of Patterson's method for the constant weight 

function. The method can be used to compute a nested sequence of quadrature formulas 

for integration with respect to any continuous probability measure on the real line with 

t*^ ■ finite moments. The advantages of the method include that it works directly with the 

moments of the underlying distribution, and that for distributions with rational moments 
the existence of the formulas can be verified by exact rational arithmetic. 
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1. Introduction 

A quadrature formula for integration with respect to the weight function p: 17 i— t- M 
takes the form 

K 

f{t)p{t)dt^Twkf{tk). (1) 

^ k=l 

^ • The weight function p is a strictly positive measurable function that is the probability 

1^ . density function of a continuous random variable with finite moments. The weights Wk 

l/^ I and nodes t^, k = 1, . . . ,K, are chosen so as to maximize the highest degree d for which the 

approximation (1) is exact for every polynomial / of degree up to d. This degree is called 
Cf^ i the degree of polynomial exactness (sometimes the degree of precision) of the formula. A 

^^ ' quadrature rule is a sequence of quadrature formulas with an increasing number of nodes 

and an increasing degree of exactness. 

For many applications it is desirable that a quadrature rule be nested, that is, that the 

node set of each formula is a subset of the node set of its successors. To obtain such a 

^ , sequence, we start with a quadrature formula with K{1) nodes, and extend it to a new 

H I formula with higher degree of exactness by adding K{2) — K{1) additional nodes. The 

" " " weights of the existing nodes may change. Repeated application of such an algorithm shall 

give rise to a nested sequence of quadrature formulas of increasing degree of polynomial 

exactness. 

In 1964, Kronrod presented a method to extend the well-known Gauss-Legendre for- 
mulas [4]. His construction adds, for any K, K + 1 nodes to the i^-node Gauss-Legendre 
formula (which has degree of exactness 2K — 1 for the constant weight function p = 1/2 
on ri = [—1,1]), extending it to a formula with 2K + 1 nodes and a degree of exact- 
ness at least 3K + 1. He also showed that this is the best possible extension (in terms 
of the achieved degree of exactness) , but he did not consider longer sequences of nested 
formulas. Patterson [5] showed that Kronrod's method can be used to obtain nested 
rules by iteratively extending the extended formulas. He also considered the constant 
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weight function over the interval [—1,1], the resulting quadrature rule is now known as 
the Gauss-Kronrod-Patterson (or GKP) rule. 

It is possible to generalize Patterson's method to obtain nested quadrature rules for 
other, non-uniform, continuous distributions with finite moments. To the best of our 
knowledge, this is very little known or used in the literature. One known example is 
the quadrature rule proposed by Genz and Keister [2] for integration with the Gaussian 
weight function. Their approach is a direct generalization of Patterson's 1964 algorithm. 
Patterson's 1989 paper [6] also considers general weight functions, but that algorithm is 
not a direct generalization of the first. 

Patterson's 1989 method requires that the distribution is given by the three-term linear 
recurrence relation satisfied by orthogonal polynomial bases with respect to the probability 
density function of the distribution. Obtaining this recurrence may itself be a difficult 
task [3], and the coefficients of the recurrence is not as commonly available information 
for known distributions as moments are. In this note we formally propose and give the 
details of an algorithm that generates nested sequences of quadrature formulas for general 
continuous distributions with finite moments. The algorithm circumvents the use of the 
three-term linear recurrence, and works with the moments of the underlying distribution 
directly. This yields a streamlined version of Patterson's algorithm, which can also be 
easily implemented. A Mathematica implementation is also provided. 

2. Extending quadrature formulas using known moments 

Our quadrature formula extension algorithm relies on the following two results. The 
first one is an immediate generalization of [4, Theorem 1]; its proof is omitted. 

Proposition 1. For every probability density function p and every set {ti, . . . ,tx} Q ^ 
of K nodes there exists unique weights wi, . . . , wk such that the quadrature formula for 
integration with respect to p has a polynomial degree of exactness at least K — 1. These 
weights are the unique solution of the linear (square) system of equations 

K 

V Wit'l = / t''p{t)dt, k = 0,...,K -1. (2) 

,^^ Jn 

Theorem 2. Let p be the probability density function of a distribution supported on Q C R 
with finite moments. Let F be a univariate polynomial of degree n with n distinct real roots, 
and suppose that there exists a polynomial G of degree p satisfying 

I F{t)G{tYp{t)dt = 0, i = 0, . . . , p - 1. (3) 

Assume further that the roots of G are all real and of multiplicity one, and distinct from 
those of F. Then there exists a quadrature formula supported on the roots of FG, whose 
degree of polynomial exactness is at least n + 2p — 1 . 

Proof. The polynomial FG has degree n + p, therefore every polynomial H of degree 
n -|- 2p — 1 can be written as F[ = QFG + R for some polynomials Q of degree p — 1 and R 
of degree n + p — 1. Consider now the quadrature formula supported on the roots of FG 
with degree of polynomial exactness at least n-\-p — l, whose existence is established in the 
previous Lemma (with n + p in place of K). Denoting by q{H) the value ^^=f WiH{ti) of 
this quadrature formula for the polynomial H we have that q{R) = J^ R{t)p{t)dt because 
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R has degree n + p — 1; q{QFG) = because q is supported on the roots of FG; and 
J^ Q{t)F{t)G{t)p{t)dt = owing to (3) and the fact that Q is of degree p — I. Hence, 

/ H{t)p{t)dt = I R{t)p{t)dt + / Q{t)F{t)G{t)p{t)dt = 
Jn Jn Jn 

= q{R) + = q{R) + q{QFG) = q{H). 

Therefore, the formula gives the correct value of J^ H{t)p{t)dt for every H of degree 
n + 2p - 1. D 

These assertions suggest the following algorithm: 

(1) Consider an arbitrary quadrature formula on n nodes, and the polynomial F whose 
roots (of multiplicity one) are the nodes of the formula. 

(2) Choose an integer p > 1, and find a degree p polynomial G that is not identically 
zero and that solves the system of equations (3). Eq. (3) is a homogeneous system 
of linear equations with one fewer variables than equations, but we can assume 
that the leading coefficient of G is 1, and turn (3) into an inhomogeneous square 
system of linear equations. 

(3) Determine the roots of G; these are the new nodes of the extended formula. Now 
find the weights of the formula by solving (2). 

By Theorem 2, the resulting formula adds new nodes to the initial formula, potentially 
increasing its degree of polynomial exactness. If p > n/2, the second formula necessarily 
has a higher degree of exactness, since an n-node formula cannot have a degree of exactness 
higher than 2n — 1. Repeated application of this algorithm may yield a nested sequence 
of quadrature formulas of increasing degree of exactness. 

As with Patterson's 1964 algorithm, this algorithm may also fail to yield a new formula 
if either the linear system (3) has no nonzero solution, or G has complex roots or real 
roots of multiplicity higher than one, or if F and G have common roots. If the algorithm 
fails for a given p, different values of p may be tried sequentially. Finally, different initial 
formulas may give rise to different sequences. 

3. Further remarks 

The original GKP formulas can be obtained using the above algorithm with Vt = [—1,1] 
and p{t) = 1/2, starting with the trivial 1-node formula (with the node at 0, with weight 
1), and taking p = n + 1 in each iteration. The number of nodes is thus 2*+^ — 1 after i 
iterations. It is not known whether this process can be repeated indefinitely, but formulas 
up to 511 nodes have been determined. The Genz-Keister formulas mentioned above can 
be obtained by successively applying the above algorithm to the normal distribution. 

It should be noted that although the computation of the nodes and weights may be 
numerically challenging (especially since the matrix of the linear system (3) is an ill- 
conditioned Hankel matrix), the existence of the solution for a given p can be decided 
using exact rational arithmetic for every distribution whose moments are rational numbers. 
In this case, the coefficients of G (if G exists) are rational numbers, and the number of 
those roots of G that are distinct from the roots of F can be determined without the 
computation of the roots of either F or G [1, Chap. 2]. 

A Mathematica implementation of the algorithm, which runs with exact rational arith- 
metic or extended precision arithmetic based on the input, is available from the authors 
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at http://users.ienis.northwestern.edu/~dpapp/, and given in the Appendix, where 
a numerical example is also provided. 

Appendix A. Mathematica code 

A simple Mathematica implementation of the quadrature rule generation algorithm con- 
sists of two functions. The first one, FormulaExtension[F,p, moments, {var, a, b},prec] 
returns the polynomial G satisfying the conditions of Theorem 2, or if such a polynomial 
does not exist. 

FormulaExtension[F_, p_, moments_, {var_, a_, b_}, prec_: \ [Infinity]] := 
Module[{G = Sum[g[i] var"i, {i, 0, p}] , 

mts = Take [moments, Exponent [F, var] + 2 p + 1], 
integrals, M, Is, success = False, ep}. 
If [prec < \[Infinity], mts = SetPrecision[mts, prec]]; 
integrals = Table [CoefficientList [F var~ipj , var]. 

Table [mts [[i+1]] , { i,0, Exponent [F var~ipj , var]}], {ipj , 0, 2 p}] ; 
M = HankelMatrix [Take [integrals , (Length [integrals] +1) /2] , 

Take [integrals , - (Length [integrals] +1) /2] ] ; 
Quiet[ls = LinearSolve [Join [Drop [M, -1], {Append [Array [0 &, p] , 1]}], 
Append [Array [0 &, p] , 1]]; 
If [Head [Is] =!= LinearSolve, 
ep = 1 s. var" Range [0, p] ; 
success = CountRoots [ep,{var,a,b}]==p && 

Resultant [F,ep, var] !=0 && Discriminant [ep, var] !=0 
] 
]; 

If [success, ep, 0] 
] 

The input arguments F and p of FormulaExtension are as in Theorem 2; moments is 
the list of known moments; the list {var,a,b} contains the symbol for the variable of 
F and the the boundaries of the interval O = [a, 6]; and prec is an optional argument 
that controls the precision of the numerical calculations. Leaving it at its default value, 
infinity, the calculations are carried out with the lowest precision of all the inputs. If 
the coefficients of F, the moments, and the endpoints of the interval [a, h] are all given 
as (symbolic) rational numbers, G all computations are carried out using exact rational 
arithmetic. 

The moments can be obtained, for instance, by symbolic integration using Mathematica. 
The moments of well-known distributions are also well-known, their precomputed values 
can be accessed through the function Moment. 

The second function needed to compute the formulas is NodesAndWeights [F, moments , var , prec] . 
which returns the nodes and weights of the formula supported on the roots of F by solving 
Eq. (2). 

NodesAndWeights [F_, moments_, var_, prec_: $MachinePrecision] : = 
Module [{w, expF = Exponent [F, var], 

nodes = var /. {ToRules [NRoots [F==0,var,PrecisionGoal->prec]] }]-, 
Transpose@{ 
nodes , 
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Array [w, expF] /. FirstSSolveSTable [ 

Sum[w[k] If [nodes [[k]]==0 && i==0, 1, nodes [[k]] "i] , {k, expF>] == 
moments [ [i+1] ] , {i,0 expF-1}] 
> 
] 

Appendix B. Example 

We illustrate the approach by computing nested quadrature formulas for the Beta(l/2,l/2) 
distribution. First, we compute the moments of the distribution: 

moms = Moment [BetaDistribution [1/2, 1/2], #] & /@ Range [0, 100]; 
Then, starting with the "empty" formula (corresponding to the polynomial constant 1), 
we find the polynomials corresponding to the successive extensions of formulas: 

F = 1; 

F = F FormulaExtension[l, 1, moms, {t, 0, 1}] 

F = F FormulaExtension[F, 2, moms, {t, 0, 1}] 

F = F FormulaExtension[F, 4, moms, {t, 0, 1}] 

F = F FormulaExtension[F, 6, moms, {t, 0, 1}] 

F = F FormulaExtension[F, 12, moms, {t, 0, 1}] 
After the last step we obtain the polynomial 

/ 1\ f y I \ f 4 ■, 19*2 3t\ f f. . 27t4 7t^ 105*2 9t 1 \ 

F{t) = it- -] [t^ -t+ — ] {t^-2f' + [t^ -3t^ + + + ■ 

V 2/ V 16/ V 16 16/ V 8 4 256 256 2048/ 

63*1" 95*9 2907*** 459*^ 1547**^ 429*^ 19305*" 1001*^ 429*^ 9* 1 

+ -^^ :::r- + ■ 



4 4 128 32 256 256 65536 32768 262144 262144 8388608 

Finally, we can obtain the formulas by running NodesAndWeights [F,moms,t ,prec] 

after each call of FormulaExtension above, with our choice of precision prec. Table 1 

shows the nodes of the obtained formulas. 

The running time of the five calls to FormulaExtension was under 0.1 seconds on a 

standard laptop with a 2.83 Ghz CPU running Mathematica 8.0.1. The CPU time required 

for the computing the weights depends on the precision. With pre c= 100 the weights of 

the five formulas were obtained in 0.4 seconds. 
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node index 

0. 3 

0.0042775693130947944277212365357185643611308627759489 5 

0.017037086855465856625128400135551316183047580495798 4 

0.038060233744356621935908405301605856588791687068179 5 

0.066987298107780676618138414623531908264298686547405 2 

0.10332332985438241771011151924935036168566203947404 5 

0. 146446609406726237799577818947575480357582031 15576 4 

0.19561928549563968029195122855091799774180314401876 5 

0.25 3 

0.30865828381745511413577000798480056661932771875719 5 

0.37059047744873961882555058118797583582546554934003 4 

0.43473690388997420422579688605225549490312964759413 5 

0.5 1 

0.56526309611002579577420311394774450509687035240587 5 

0.62940952255126038117444941881202416417453445065997 4 

0.69134171618254488586422999201519943338067228124281 5 

0.75 3 

0.80438071450436031970804877144908200225819685598124 5 

0.85355339059327376220042218105242451964241796884424 4 

0.89667667014561758228988848075064963831433796052596 5 

0.93301270189221932338186158537646809173570131345260 2 

0.96193976625564337806409159469839414341120831293182 5 

0.98296291314453414337487159986444868381695241950420 4 

0.99572243068690520557227876346428143563886913722405 5 

1. 3 

Table 1. Nodes of a nested quadrature rule for the Beta(l/2,l/2) distri- 
bution. The index column shows the smallest k for which the node appears 
in the A;th formula. The weights of the formulas are omitted for brevity. 



